% Migration Data analysis
clc
clear

Y = cellstr(char('1998','1999','2000','2001','2002','2003','2004','2005','2006','2007','2008','2009'));
R = cellstr(char('1.5','2','2.5','3','3.5','4','4.5'));
SG = cellstr(char('Sample 1: [$50k,+)','Sample 2A: [$50k, 250k)','Sample 2B: [$250k, 1000k)','Sample 2C: [$1000k, +)'));
Markers=['o','+','s','d','*','x','v','^','<','>','p','h','.',...
'+','*','o','x','^','<','h','.','>','p','s','d','v',...
'o','x','+','*','s','d','v','^','<','>','p','h','.'];
linestyles = cellstr(char('-',':','-.','--','-',':','-.','--','-',':','-',':',...
'-.','--','-',':','-.','--','-',':','-.'));

load duration1yrs_TP06YRJfinR5G50MR7.mat durationCount NUMBER_OF_YEAR RATE_LIST_LENGTH ;
durationCountSample = zeros(length(SG),NUMBER_OF_YEAR,RATE_LIST_LENGTH-1,RATE_LIST_LENGTH );
durationCountSample(1,:,:,:) = durationCount;
clear durationCount

load duration1yrs_TP06YRJfinR5G50250MR7.mat durationCount
durationCountSample(2,:,:,:) =durationCount;
clear durationCount

load duration1yrs_TP06YRJfinR5G2501000kMR7.mat durationCount
durationCountSample(3,:,:,:) =durationCount;
clear durationCount

load duration1yrs_TP06YRJfinR5G1000kMR7.mat durationCount 
durationCountSample(4,:,:,:) =durationCount;
clear durationCount




%% Unconditional Transition Matrix


stdProbY = zeros (length(SG), NUMBER_OF_YEAR, RATE_LIST_LENGTH-1, RATE_LIST_LENGTH);
stdProbAvgY = zeros (length(SG), RATE_LIST_LENGTH-1, RATE_LIST_LENGTH);


for s = 1:length(SG)
    
    xx = zeros (RATE_LIST_LENGTH-1, RATE_LIST_LENGTH);
    for i = 1: NUMBER_OF_YEAR
        x(:, :) = durationCountSample(s,i,:,:);
        xx = xx + x;
        for d = 1: RATE_LIST_LENGTH-1
            stdProbY (s, i, d, :) = x(d,:)./ sum(x(d,:),2);
        end
    end
    
    for d = 1: RATE_LIST_LENGTH-1
        stdProbAvgY (s, d, :) = xx(d,:)./ sum(xx(d,:),2);
    end
    
    clear  xx
end


% Draw the figure to check the downgrade issue for rating 4
% Will draw figure of migration of rating 1.5-3 along the time.

f = zeros(NUMBER_OF_YEAR,RATE_LIST_LENGTH);
Linecolors=jet(length(Y));
for s = 1:length(SG)
    figure
    text(0.5,0.5,'SG(s'); % it doesn't work.
    for j = 1:4 % for each rating 1.5 -3
        subplot(2,2,j)
        for i = 1: NUMBER_OF_YEAR
            f(i,:) = stdProbY(s,i, j, :);
            hold on
            plot(f(i,j+1:RATE_LIST_LENGTH)',[linestyles{1} Markers(i)],'LineWidth',2,'Color',Linecolors(i,:),'MarkerSize',10);
        end
        hold off
        legend(Y)
        set(gca,'XTick',1:length(R))
        set(gca,'XTickLabel',[R(j+1:length(R)); 'Default']);
        xlabel('Ratings');
        title(strcat('Downgrade Migration from Rating -', R(j)));
        grid on
    end
    
end

% % To check only column of rating 4 probability. Not good effect.
% ff = zeros(NUMBER_OF_YEAR,strmatch({'3.5'},R));
% Linecolors=jet(length(Y));
% figure
% for i = 1: NUMBER_OF_YEAR
% ff(i,:) = stdProbY(i, 1:strmatch({'3.5'},R), strmatch({'4.5'},R)-1); % take column rating 4 for first 5 rows (rating 1.5-3.5)
% hold on
% plot(ff(i,:),[linestyles{1} Markers(i)],'LineWidth',2,'Color',Linecolors(i,:),'MarkerSize',10);
% end
% hold off
% legend(Y)
% set(gca,'XTick',1:strmatch({'3.5'},R))
% set(gca,'XTickLabel',R(1:strmatch({'3.5'},R)));
% xlabel('Ratings');
% ylabel('Downgrade Probability')
% grid on

%% Diagnoal

diagProbY = zeros (length(SG),NUMBER_OF_YEAR, RATE_LIST_LENGTH-1);
Linecolors=jet(length(R));
figure;

for s = 1: length(SG)
   
    
    for d = 1: RATE_LIST_LENGTH-1
        for q = 1 : NUMBER_OF_YEAR
            diagProbY (q,d) = stdProbY (s, q, d, d);
        end
        subplot(2,2,s);
        hold on
        plot(diagProbY(:,d),[linestyles{1} Markers(d)],'LineWidth',2,'Color',Linecolors(d,:),'MarkerSize',10);
    end
    hold off
    
    set(gca,'XLim',[1 length(diagProbY)]);
    set(gca,'XTickLabel',Y);
    xlabel('Years');
    ylabel('Non-Transition Frequency');
    legend(R,'Location','NorthWest');
    title(SG(s));
    grid on
end

%% Coefficient of Variation
% coefficient of variation is defined as the ratio of the standard deviation divided by the absolute vaule of mean.

CV = zeros (length(SG), RATE_LIST_LENGTH-1, RATE_LIST_LENGTH);
AvgCV = zeros(length(SG),1);
filename='D:\My Documents\CIRANO\Writing\Ch1_Data Statistics.2011nov9\tbldata.xlsx';
sheetlist=strcat('coef var');

for s = 1: length(SG)
    
    xx = zeros (RATE_LIST_LENGTH-1, RATE_LIST_LENGTH);
    miu(:,:) = stdProbAvgY (s, :, :);
    
    for i=1:NUMBER_OF_YEAR
        x(:,:) = stdProbY (s, i, :, :);
        xx(:,:)= xx(:,:) +((x-miu).^2); %.^2 is to square elements of matrix.
    end
    
    sigma(:,:)=(xx./12).^(0.5);
    xxx =sigma./abs(miu(:,:));
    CV(s,:,:)= xxx;
    AvgCV(s,:) = mean(mean(xxx));
    y = 1+(s-1)*7;
    xlswrite(filename, xxx, sheetlist, strcat('A',num2str(s+(s-1)*7)));
    
    clear x xx xxx
    
end

figure
for s = 1: length(SG)
    subplot(2,2,s);
    xxx (:,:)= CV(s,:,:); 
    bar3(xxx',0.35,'detached');   
    set(gca,'XDir','reverse');
    set(gca, 'XTickLabel', R);
    set(gca, 'YTickLabel', [R;'D']);
    set(gca, 'Ylim',[0.5, 8.5]);
    set(gca,'ZLim',[0, 0.8]);
    xlabel('Initial Rating');
    ylabel('End Rating');
    zlabel('Coeffient Variation');
    title (SG(s))
end

%% Downgrade and upgrade along the time
updownCount = zeros(length(SG),length(Y),2); % first column is up, and 2nd column is down
updownPcg = zeros(length(SG),length(Y),2); % first column is up, and 2nd column is down

for s = 1: length(SG)
    for i = 1: length(Y)
        xx(:,:) = durationCountSample(s,i,:,:);
        
        % count downgrade number
        for r = 1: length(R)
            for c = r+1: length(R)+1
                updownCount (s,i,2) = updownCount(s,i,2) + xx(r,c);   
            end
        end
        
        % count upgrade number
        for r = 2: length(R)
            for c = 1: r-1
                updownCount (s,i,1)= updownCount(s,i,1) + xx(r,c);                
            end
        end
        
        % count percentage
        updownPcg(s,i,1) = updownCount(s,i,1)/ sum(sum(xx));
        updownPcg(s,i,2) = updownCount(s,i,2)/ sum(sum(xx));
        
    end
end

figure
for s = 1: length(SG)
          
    subplot(2,2,s);
    y1(:,:) = updownCount(s,:,:);
    y2(:,:) = updownPcg(s,:,:);
    
    x = 1998:1:2009;
    [AX,H1,H2]=plotyy(x,y1,x,y2,'bar','plot');
    colormap(summer);
    set(AX,'xlim',[1997,2010]);
    xlabel('Years');
    set(get(AX(1),'Ylabel'),'String','Number') 
    set(get(AX(2),'Ylabel'),'String','Percentage') 
    set(AX(1),'YAxisLocation','left','YTickMode','auto')
    set(AX(2),'YAxisLocation','right','YTickMode','auto')
    set(H2(1), 'color','r','marker', 'o','Linewidth',2,'MarkerSize',10);
    set(H2(2), 'color','b','marker', 's','Linewidth',2,'MarkerSize',10);     
    set(gca,'XTickLabel',Y);
    axes(AX(1)),
    legend('Upgrade','Downgrade',2)
    axes(AX(2)),
    legend('Upgrade %','Downgrade %')      
    title(SG(s));
    grid on
    
end



% % evolution of rating universe over time-----------------------------------
% 
% numIssuers = zeros (NUMBER_OF_YEAR * Quarter_LIST_LENGTH, 1);
% 
% for i = 1: NUMBER_OF_YEAR
%     for j = 1: Quarter_LIST_LENGTH
%         x(:, :) = durationCount(i,j,:,:);
%         numIssuers(4*(i-1)+j,:) = sum(sum(x));
%     end
% end
% numIssuersTS = timeseries(numIssuers(:,:), 1: length(numIssuers), 'name', 'Issuers Quarter');
% 
% figure
% plot (numIssuersTS,'-*');
% set(gca,'XTick',1:4:length(numIssuers));
% labels = quaterlabels(1998, length(numIssuers));
% labels = labels(1:4:length(numIssuers));
% set(gca,'XTickLabel',labels);
% %grid minor;
% xlabel('Quarter');
% 
% % Average Rating Distribution----------------------------------------------
% % Evolution of Rating Over Time
% 
% rateDist = zeros (RATE_LIST_LENGTH-1, 1);
% rateDist06Before = zeros (RATE_LIST_LENGTH-1, 1);
% rateOverTime = zeros( NUMBER_OF_YEAR * Quarter_LIST_LENGTH , RATE_LIST_LENGTH-1 );
% for i = 1: NUMBER_OF_YEAR
%     for j = 1: Quarter_LIST_LENGTH
%         x(:, :) = durationCount(i,j,:,:);
%         rateDist = sum(x,2)+ rateDist;  % exclude withdrawn, how many obligors quarters totally.  
%         if i==8 
%             rateDist06Before = rateDist; % exclude withdrawn, how many obligors quarters totally before 2006q1.
%         end          
%         rateOverTime (4*(i-1)+j,:)= sum(x,2)';
%     end
% end
% 
% rateDist06After = rateDist - rateDist06Before; % exclude withdrawn, how many obligors quarters totally after 2006q1.
% 
% rateDistAvg = rateDist./(NUMBER_OF_YEAR*Quarter_LIST_LENGTH);
% rateDistAvgBefore2006 = rateDist06Before./(8*Quarter_LIST_LENGTH);
% rateDistAvgAfter2006 = rateDist06After./(4*Quarter_LIST_LENGTH);
% 
% 
% figure;
% bar (rateDistAvg);
% set(gca,'XTickLabel',{'1.5','2','2.5','3','3.5','4','4.5' });
% xlabel('Ratings');
% ylabel('Obligors Quarter');
% title('Average Rating Distribution');
% 
% figure;
% bar (rateDistAvgBefore2006);
% set(gca,'XTickLabel',{'1.5','2','2.5','3','3.5','4','4.5' });
% xlabel('Ratings');
% ylabel('Obligors Quarter');
% title('Average Rating Distribution Before 2006Q1');
% 
% figure;
% bar (rateDistAvgAfter2006);
% set(gca,'XTickLabel',{'1.5','2','2.5','3','3.5','4','4.5' });
% xlabel('Ratings');
% ylabel('Obligors Quarter');
% title('Average Rating Distribution After 2006Q1');
% 
% figure;
% plot (rateOverTime);
% set(gca,'XTick',1:4:length(numIssuers));
% set(gca,'XTickLabel',labels);
% xlabel('Quarter');
% ylabel('Obligors Quarter');
% legend('1.5','2','2.5','3','3.5','4','4.5','Location','NorthWest');
% 
% 
% % Default Freqency Evolution-----------------------------------------------
% % Evolution of Default Frequency Over Time and Ratings
% 
% defaultFreq = zeros (NUMBER_OF_YEAR * Quarter_LIST_LENGTH, 1);
% defaultFreqOverRate = zeros (NUMBER_OF_YEAR * Quarter_LIST_LENGTH, RATE_LIST_LENGTH-1);
% 
% for i = 1: NUMBER_OF_YEAR
%     for j = 1: Quarter_LIST_LENGTH
%         x(:, :) = durationCount(i,j,:,:);
%         defaultFreq(4*(i-1)+j,:) = sum(x(:,RATE_LIST_LENGTH)) / sum(sum(x));
%         defaultFreqOverRate(4*(i-1)+j,:) = ( x(:, RATE_LIST_LENGTH) ./ sum(x, 2) )';
%     end
% end
% 
% defaultFreqTS = timeseries(defaultFreq(:,:), 1: length(defaultFreq), 'name', 'Default Freqency');
% 
% figure;
% plot (defaultFreqTS,'-*');
% set(gca,'XTick',1:4:length(defaultFreq));
% set(gca,'XTickLabel',labels);
% grid minor;
% xlabel('Quarter');
% 
% figure;
% plot (defaultFreqOverRate);
% set(gca,'XTick',1:4:length(defaultFreq));
% set(gca,'XTickLabel',labels);
% xlabel('Quarter');
% ylabel('Default Freqency');
% legend('1.5','2','2.5','3','3.5','4','4.5','Location','NorthWest');




% % Metrics Calculation
% % SVD (define singular value of matrix as mobility norm, here we use average of singular
% % value instead of largest one which is argued in Jafry and Schuerman)
% % D3 (risk-adjusted indices)
% 
% % NOTICE !! Here all the metrics is distance to average matrix (all sample)
% mobilityNorm = zeros (NUMBER_OF_YEAR * Quarter_LIST_LENGTH,1);
% mobilityNormRate = zeros (NUMBER_OF_YEAR * Quarter_LIST_LENGTH,1);
% SVD = zeros(NUMBER_OF_YEAR * Quarter_LIST_LENGTH,1);
% SVDRate = zeros(NUMBER_OF_YEAR * Quarter_LIST_LENGTH,1);
% 
% L1 = zeros (NUMBER_OF_YEAR * Quarter_LIST_LENGTH,1);
% L2 = zeros (NUMBER_OF_YEAR * Quarter_LIST_LENGTH,1);
% WAD = zeros (NUMBER_OF_YEAR * Quarter_LIST_LENGTH,1);
% NAD = zeros (NUMBER_OF_YEAR * Quarter_LIST_LENGTH,1);
% D3 = zeros (NUMBER_OF_YEAR * Quarter_LIST_LENGTH,1);
% D1 = zeros (NUMBER_OF_YEAR * Quarter_LIST_LENGTH,1);
% D2 = zeros (NUMBER_OF_YEAR * Quarter_LIST_LENGTH,1);
% D4 = zeros (NUMBER_OF_YEAR * Quarter_LIST_LENGTH,1);
% D1sqr = zeros (NUMBER_OF_YEAR * Quarter_LIST_LENGTH,1);
% D2sqr = zeros (NUMBER_OF_YEAR * Quarter_LIST_LENGTH,1);
% D3sqr = zeros (NUMBER_OF_YEAR * Quarter_LIST_LENGTH,1);
% D4sqr = zeros (NUMBER_OF_YEAR * Quarter_LIST_LENGTH,1);
% 
% 
% rateMatrix = stdProbQ (:, :, 1:RATE_LIST_LENGTH-1);
% rateMatrixAvg = stdProbAvgQ (:, 1:RATE_LIST_LENGTH-1);
% a = mean(svd(stdProbAvgQ));
% aRate = mean(svd(rateMatrixAvg));
% 
% for i = 1: size(stdProbQ, 1)
%     aa(:,:) = stdProbQ (i, :, :);
%     aaRate (:,:) = rateMatrix (i, :, :);
%     mobilityNorm (i,1) = mean(svd(aa));
%     mobilityNormRate(i,1) = mean(svd(aaRate));
%     SVD(i,1) = mobilityNorm (i,1) - a; %SVD
%     SVDRate(i,1) = mobilityNormRate (i,1) - aRate;
%     
%     for j = 1 : RATE_LIST_LENGTH-1
%         for k = 1 : RATE_LIST_LENGTH
%             
%             L1(i) = L1(i) + abs(aa(j,k) - stdProbAvgQ(j,k)); % L1
%             L2(i) = L2(i) + (aa(j, k) - stdProbAvgQ(j, k))^2;
%             WAD(i)= WAD(i) +  abs(aa(j,k) - stdProbAvgQ(j,k))* aa(j,k); %WAD
%             
%             
%             if aa(j,k)~=0
%                 NAD(i) =  NAD (i) + abs(aa(j,k) - stdProbAvgQ(j,k))/aa(j,k); %NAD
%                 mD2 = (j - k)* (aa(j, k) - stdProbAvgQ(j, k)) / aa(j, k);
%                 mD4 = (j - k)* sign(aa(j, k) - stdProbAvgQ(j, k))*(aa(j, k) - stdProbAvgQ(j, k))^2 / aa(j, k);
%             end
%             
%             
%             
%             mD3 = (j - k)* sign(aa(j, k) - stdProbAvgQ(j, k))*(aa(j, k) - stdProbAvgQ(j, k))^2;
%             mD1 = (j - k)* (aa(j, k) - stdProbAvgQ(j, k));
%             
%             if k~= RATE_LIST_LENGTH %D3
%                 D3(i) = D3(i) + mD3;
%                 D1(i) = D1(i) + mD1;
%                 D2(i) = D2(i) + mD2;
%                 D4(i) = D4(i) + mD4;
%                 D1sqr(i)= D1sqr(i)+ mD1;
%                 D2sqr(i)= D1sqr(i)+ mD2;
%                 D3sqr(i)= D1sqr(i)+ mD3;
%                 D4sqr(i)= D1sqr(i)+ mD4;
%                 
%             else
%                 D3(i) = D3(i) + mD3*RATE_LIST_LENGTH;
%                 D1(i) = D1(i) + mD1*RATE_LIST_LENGTH;
%                 D2(i) = D2(i) + mD2*RATE_LIST_LENGTH;
%                 D4(i) = D4(i) + mD4*RATE_LIST_LENGTH;                
%                 D1sqr(i) = D1sqr(i)+ mD1*RATE_LIST_LENGTH^2;
%                 D2sqr(i) = D2sqr(i)+ mD2*RATE_LIST_LENGTH^2;
%                 D3sqr(i) = D3sqr(i)+ mD3*RATE_LIST_LENGTH^2;
%                 D4sqr(i) = D4sqr(i)+ mD4*RATE_LIST_LENGTH^2;
%             end
%             
%             
%         end
%     end
%     L2(i) = sqrt(L2(i)); %L2
% end
% 
% 
% % 
% % figure;
% % plot (mobilityNorm);
% % set(gca,'XTick',1:4:length(mobilityNorm));
% % set(gca,'XTickLabel',labels);
% % xlabel('Quarter');
% % ylabel('Mobility Norm SVD');
% 
% figure;
% subplot(4,2,1); plot (SVD);
% set(gca,'XTick',1:4:length(mobilityNorm));
% set(gca,'XTickLabel',labels);
% xlabel('Quarter');
% ylabel('SVD');
% 
% subplot(4,2,2);plot (D3);
% set(gca,'XTick',1:4:length(D3));
% set(gca,'XTickLabel',labels);
% xlabel('Quarter');
% ylabel('Risk-adjusted Metric: D');
% 
% subplot(4,2,3);plot (L1);
% set(gca,'XTick',1:4:length(L1));
% set(gca,'XTickLabel',labels);
% xlabel('Quarter');
% ylabel('L1');
% 
% subplot(4,2,4);plot (L2);
% set(gca,'XTick',1:4:length(L2));
% set(gca,'XTickLabel',labels);
% xlabel('Quarter');
% ylabel('L2');
% 
% subplot(4,2,5);plot (WAD);
% set(gca,'XTick',1:4:length(WAD));
% set(gca,'XTickLabel',labels);
% xlabel('Quarter');
% ylabel('WAD');
% 
% subplot(4,2,6);plot (NAD);
% set(gca,'XTick',1:4:length(NAD));
% set(gca,'XTickLabel',labels);
% xlabel('Quarter');
% ylabel('NAD');
% 
% subplot(4,2,7);plot (D1);
% set(gca,'XTick',1:4:length(WAD));
% set(gca,'XTickLabel',labels);
% xlabel('Quarter');
% ylabel('D1');
% 
% subplot(4,2,8);plot (D1sqr);
% set(gca,'XTick',1:4:length(NAD));
% set(gca,'XTickLabel',labels);
% xlabel('Quarter');
% ylabel('D1sqr');




% 

% % Distribution of PD and inverse normal CDF of PD--------------------------
% % Distribution of SVD and stdSVD
% % Distribution of D3
% % Distribution of L1
% % Distribution of WAD, NAD
% % 
% 
% % Compare the actual distribution with normal distribution.
% 
% defFreqInverse = norminv(defaultFreq);
% stdSVD = -zscore(SVD);
% 
% 
% figure
% subplot (1, 2, 1);
% probplot(defaultFreq);
% legend('Normal','PD','Location','NW')
% subplot (1, 2, 2);
% probplot(defFreqInverse);
% legend('Normal','Inverse Normal CDF of PD','Location','NW')
% 
% figure
% 
% subplot (4, 2, 1);
% probplot(L1);
% legend('Normal','L1','Location','NW');
% subplot (4, 2, 2);
% probplot(L2);
% legend('Normal','L2','Location','NW')
% subplot (4, 2, 3);
% probplot(WAD);
% legend('Normal','WAD','Location','NW')
% subplot (4, 2, 4);
% probplot(NAD);
% legend('Normal','NAD','Location','NW')
% subplot (4, 2, 5);
% probplot(SVD);
% legend('Normal','SVD','Location','NW');
% subplot (4, 2, 6);
% probplot(D3);
% legend('Normal','D3','Location','NW')
% subplot (4, 2, 7);
% probplot(D1);
% legend('Normal','D1','Location','NW');
% subplot (4, 2, 8);
% probplot(D1sqr);
% legend('Normal','D1sqr','Location','NW')
% 
% % probplot(defaultFreq)
% % p = mle(defaultFreq,'dist','tlo');
% % t = @(defaultFreq,mu,sig,df)cdf('tlocationscale',defaultFreq,mu,sig,df);
% % h = probplot(gca,t,p);
% % set(h,'color','r','linestyle','-')
% 
% 
% % chi-square goodness-of-fit test for testing normal distribution of PD and
% % inverse normal CDF of PD
% [h(1),p(1)] = chi2gof(defaultFreq);
% [h(2),p(2)] = chi2gof(defFreqInverse);
% [h(3),p(3)] = chi2gof(L1);
% [h(4),p(4)] = chi2gof(L2);
% [h(5),p(5)] = chi2gof(WAD);
% [h(6),p(6)] = chi2gof(NAD);
% [h(7),p(7)] = chi2gof(SVD);
% [h(8),p(8)] = chi2gof(D3);
% [h(9),p(9)] = chi2gof(D1);
% [h(10),p(10)] = chi2gof(D1sqr);
% 
% 
% % draw graph for SVD, D1-D4 , D1_n^2 - D4_n^2, (L1,WAD,NAD) distance to average TM
% t = 1:1:NUMBER_OF_YEAR * Quarter_LIST_LENGTH;
% 
% figure;
% plot(t,0); hold on % add recession figure on previous graph
% ha = area([41 47], [0.03 0.03]);
% set(ha,'BaseValue',-0.05);
% set(ha,'FaceColor',[0.83 0.82 0.78]);
% set(ha,'LineStyle','none')
% grid on
% set(gca,'Layer','top'); hold on % add metrics on previous graph
% hb = plot (SVD,'-+');
% set(gca,'XTick',1:4:length(SVD));
% set(gca,'XTickLabel',labels);
% xlabel('Quarter');
% ylabel('Mobility Norm SVD');
% grid on
% 
% figure;
% plot(t,0); hold on % add recession figure on previous graph
% ha = area([41 47], [5 5]);
% set(ha,'BaseValue',-4);
% set(ha,'FaceColor',[0.83 0.82 0.78]);
% set(ha,'LineStyle','none')
% grid on
% set(gca,'Layer','top'); hold on % add metrics on previous graph
% hb = plot(t,zscore(D1),'-+',t,zscore(D2),'-o',t,zscore(D3),'-*',t,zscore(D4),'-x');
% legend(hb,'D1','D2','D3','D4','Location','NorthWest');
% set(gca,'XTick',1:4:length(numIssuers));
% set(gca,'XTickLabel',labels);
% xlabel('Quarter');
% ylabel('Risk-adjusted Difference Metrics D1-D4');
% grid on
% 
% 
% figure;
% plot(t,0); hold on % add figure on previous graph
% ha = area([41 47], [4 4]);
% set(ha,'BaseValue',-3);
% set(ha,'FaceColor',[0.83 0.82 0.78]);
% set(ha,'LineStyle','none')
% grid on
% set(gca,'Layer','top'); hold on
% hb = plot(t,zscore(D1sqr),'-+',t,zscore(D2sqr),'-o',t,zscore(D3sqr),'-*',t,zscore(D4sqr),'-x');
% legend(hb,'D1(n^2)','D2(n^2)','D3(n^2)','D4(n^2)','Location','NorthWest');
% set(gca,'XTick',1:4:length(numIssuers));
% set(gca,'XTickLabel',labels);
% xlabel('Quarter');
% ylabel('Risk-adjusted Difference Metrics D1_n^2-D4_n^2');
% grid on
% 
% figure;
% plot(t,0); hold on % add figure on previous graph
% ha = area([41 47], [4 4]);
% set(ha,'BaseValue',-2);
% set(ha,'FaceColor',[0.83 0.82 0.78]);
% set(ha,'LineStyle','none')
% grid on
% set(gca,'Layer','top'); hold on
% hb = plot(t,zscore(L1),'-+',t,zscore(WAD),'-o',t,zscore(NAD),'-*');
% legend(hb,'L1','WAD','NAD','Location','NorthWest');
% set(gca,'XTick',1:4:length(numIssuers));
% set(gca,'XTickLabel',labels);
% xlabel('Quarter');
% ylabel('Risk-adjusted Difference Metrics L1,WAD,NAD');
% grid on
% 
% 
% % draw the graph of cell prob along time
% % We did the same cell as paper Xing (2010)
% % ------------------
% % 1 AAA 1.5
% % 2 AA  2
% % 3 A   2.5
% % 4 BBB 3
% % 5 BB  3.5
% % 6 B   4
% % 7 C   4.5
% % 8 D   D
% %-------------------
% % the cell to check include: 1-1, 5-5, 1-8, 5-8, 2-6, 3-4
% 
% figure;
% subplot(3,2,5); plot (stdProbQ(:,2,2));
% legend('Rating 2 to 2 ');
% title('Observed Transition Prob along Time',... 
%   'FontWeight','bold')
% set(gca,'XTick',1:4:length(numIssuers));
% labels = quaterlabels(1998, length(numIssuers));
% labels = labels(1:4:length(numIssuers));
% set(gca,'XTickLabel',labels);
% xlabel('Quarter');
% grid on













% graph for regression

% SVD regression


% Analyze the trend situation.
% figure
% detrend_dp=detrend(defaultFreq);
% trend = defaultFreq - detrend_dp;
% plot(defaultFreq);
% legend('Original Data','Location','northwest');
% hold on
% plot(trend,':r')
% plot(detrend_dp,'m')
% plot(zeros(size(defaultFreq)),':k')
% legend('Original Data','Trend','Detrended Data',...
%        'mean(Detrended)','Location','northwest')
% xlabel('Quarters'); 
% ylabel('DP');



% Test against the standard normal:chi2gof(SVD,'cdf',@normcdf)
% Use lillietest to determine if SVD follows a normal distribution: lillietest(SVD)
% The Lilliefors test is a 2-sided goodness-of-fit test suitable when a fully-specified null distribution is unknown and its parameters must be estimated. This is in contrast to the one-sample Kolmogorov-Smirnov test (see kstest), which requires that the null distribution be completely specified.





